C ***********************************************************************
C *            Subroutine LEBEM3 by Stephen Kirkup                      *
C ***********************************************************************
C
C Website : www.boundary-element-method.com (LAPLACE/BEMLAP)
C
C This subroutine computes the solution to the three-dimensional 
C Laplace equation
C                  __ 2                
C                  \/    {\phi}   =  0   
C
C in the domain exterior to a closed boundary.
C
C The boundary (S) is defined (approximated) by a set of planar 
C triangular panels. The domain of the equation is within the
C boundary.
C
C The boundary condition may be Dirichlet, Robin or Neumann. It is 
C assumed to have the following general form
C
C           {\alpha}(q) {\phi}(q) + {\beta}(q) v(q) = f(q)
C    
C where {\phi}(q) is the solution at the point q on S, v(q) is the 
C derivative of {\phi} with respect to the outward normal to S at q and
C {\alpha}, {\beta} and f are real-valued functions defined on S. The
C functions {\alpha} and {\beta} must be specified to define the nature
C of the boundary condition. Important examples are {\alpha}=1, 
C {\beta}=0 which is equivalent to a Dirichlet boundary condition and 
C {\alpha}=0, {\beta}=1 which is equivalent to a Neumann boundary 
C condition. The specification of f completes the definition of the 
C boundary condition.
C
C
C How to use the subroutine
C -------------------------
C
C The following diagram shows how the subroutine is to be used. 
C
C                                   ....................................
C                                   :                                  :
C                                   :                                  :
C      ----------------------       :     --------------------------   :
C      |                    |       :     |                        |   :
C      |   MAIN PROGRAM     |------->-----|      LEBEM3            |   :
C      | (e.g.LEBEM3_t.for) |       :     |                        |   :
C      |                    |       :     --------------------------   :
C      ----------------------       :                 |                :
C                                   :                 >                :
C                                   :                 |                :
C                                   :      ------------------------    :
C          Package ---------------->:      | subordinate routines |    :
C                                   :      ------------------------    :
C                                   :                                  :
C                                   :      (this file)                 :  
C                                   :..................................:
C                                  /         |                 |
C                               |_           >                 >
C                              /             |                 |
C             ................       ................   ................  
C             :              :       :   --------   :   :  --------    : 
C             : (geom3d.for) :---<---:   | L3LC |   :   :  | GLS |    : 
C             :              :       :   --------   :   :  --------    :  
C             :..............:       : -------------:   : -------------:  
C                                    : |subordinate|:   : |subordinate|: 
C                                    : | routines  |:   : | routines  |:  
C                                    : -------------:   : -------------: 
C                                    :              :   :              : 
C                                    : (l3lc.for)   :   : (gls.for)   :
C                                    :..............:   :..............:
C                                    
C
C The contents of the main program must be linked to LEBEM3.FOR, L3LC.FOR
C and GLS.FOR.
C
C Method of solution
C ------------------
C 
C In the main program, the boundary must be described as a set of
C  panels. The panels are defined by three indices (integers) which
C  label a node or vertex on the boundary. The data structure VERTEX 
C  lists and enumerates the coordinates of the vertices, the data 
C  structure SELV defines each panel by indicating the labels for
C  the three nodes that are its vertices and hence enumerates the
C  panels.
C The boundary solution points (the points on the boundary at which 
C  {\phi} (SPHI) and d {\phi}/dn (SVEL) are returned) are at the centres
C  of the panels. The boundary functions {\alpha} (SALPHA), {\beta} 
C  (SBETA) and f (SF) are also defined by their values at the centres
C  of the panels.
C Normally a solution in the domain is required. By listing the 
C  coordinates of all the exterior points in PENT, the subroutine
C  returns the value of {\phi} at these points in PEPHI.
C
C
C Format and parameters of the subroutine
C ---------------------------------------
C
C The subroutine has the form
C
C      SUBROUTINE LEBEM3(MAXNV,NV,VERTEX,MAXNSE,NSE,SELV,
C     *                 MAXNPE,NPE,PENT,
C     *                 SALPHA,SBETA,SF,SIPHI,SIVEL,PDEPHI,
C     *                 LSOL,LVALID,EGEOM,
C     *                 SPHI,SVEL,PEPHI,
C     *                 WKSPC1,WKSPC2,WKSPC3,WKSPC4,
C     *                 WKSPC5,WKSPC6,WKSPC7)

C The parameters to the subroutine
C ================================

C Geometry of the boundary S (input)
C integer MAXNV: The limit on the number of vertices of the polygon
C  that defines (approximates) S. MAXNV>=4.
C integer NV: The number of vertices on S. 4<=NV<=MAXNV.
C real VERTEX(MAXNV,3): The coordinates of the vertices. VERTEX(i,1),
C  VERTEX(i,2), VERTEX(i,3) are the x,y,z coordinates of the i-th 
C  vertex.
C integer MAXNSE: The limit on the number of panels describing S.
C  MAXNSE>=4.
C integer NSE: The number of panels describing S. 4<=NSE<=MAXNSE.
C integer SELV(MAXNSE,3): The indices of the three vertices defining
C  each panel. The i-th panel have vertices 
C  (VERTEX(SELV(i,1),1),VERTEX(SELV(i,1),2)),VERTEX(SELV(i,1),3)),
C  (VERTEX(SELV(i,2),1),VERTEX(SELV(i,2),2)),VERTEX(SELV(i,2),3)) and
C  (VERTEX(SELV(i,3),1),VERTEX(SELV(i,3),2)),VERTEX(SELV(i,3),3)).

C Interior points at which the solution is to be observed (input)
C integer MAXNPE: Limit on the number of points exterior to the
C  boundary. MAXNPE>=1.
C integer NPE: The number of exterior points. 0<=NPE<=MAXNPE.
C real PENT(MAXNPE,3). The coordinates of the exterior point.
C  PENT(i,1),PENT(i,2),PENT(i,3) are the x,y,z coordinates of the i-th
C  point.

C The boundary condition ({\alpha} phi + {\beta} v = f) (input)
C real SALPHA(MAXNSE): The values of {\alpha} at the centres
C  of the panels.
C real SBETA(MAXNSE): The values of {\beta} at the centres
C  of the panels.
C real SF(MAXNSE): The values of "f" at the centres of the 
C  panels.

C Incident field
C real SIPHI(MAXNSE): The incident velocity potential at the
C  centres of the panels
C real SIVEL(MAXNSE): The derivative of the incident velocity 
C  centres of the panels
C real PDEPHI(MAXNPE): The incident velocity potential at the chosen
C  exterior points

C Validation and control parameters (input) 
C logical LSOL: A switch to control whether the particular solution is
C  required.
C logical LVALID: A switch to enable the choice of checking of 
C  subroutine parameters.
C real EGEOM: The maximum absolute error in the parameters that
C  describe the geometry.

C Solution (output)
C real SPHI(MAXNSE): The velocity potential ({\phi}) at the 
C  centres of the boundary panels.
C real SVEL(MAXNSE): The velocity (v or d{\phi}/dn where n is
C  the outward normal to the boundary) at the centres of the boundary 
C  panels.
C real PEPHI(MAXNPE): The velocity potential ({\phi}) at the 
C  exterior points.

C Working space
C real    WKSPC1(MAXNSE,MAXNSE)
C real    WKSPC2(MAXNSE,MAXNSE)
C real    WKSPC3(MAXNPE,MAXNSE)
C real    WKSPC4(MAXNPE,MAXNSE)
C real    WKSPC5(MAXNSE)
C real    WKSPC6(MAXNSE)
C logical WKSPC7(MAXNSE)

C Notes on the geometric parameters
C ---------------------------------
C (1) Each of the vertices listed in VERTEX must be distinct points
C  with respect to EGEOM.
C (2) The boundary must be complete and closed. 
C (3) The indices of the nodes listed in SELV must be such that they
C  are ordered counter-clockwise around the boundary, when viewed
C  from just outside the boundary at the panel.
C (4) The largest panel must be no more than 10x the area of the
C  smallest panel.

C Notes on the exterior points 
C ----------------------------
C (1) The points in PENT should lie within the boundary, as defined
C  by the parameters VERTEX and SELV. Any point lying outside the 
C  boundary will return a corresponding value in PEPHI that is near
C  zero. This property can be if a useful check.

C Notes on the boundary condition
C -------------------------------
C (1) For each i=1..NSE, it must not be the case that both of SALPHA(i)
C  and SBETA(i) are zero

C External modules in external files
C ==================================
C subroutine L3LC: Returns the individual discrete Laplace integral
C  operators. (in file L3LC.FOR)
C subroutine GLS: Solves a general linear system of equations. 
C  (in file GLS.FOR)

C External modules provided in the package (this file)
C ====================================================
C subroutine GLT7: Returns the points and weights of the 7-point Gauss-
C  Legendre quadrature rule on the standard triangle.
C real function FNSQRT(X): real X : Returns the square root of X.
C real function FNEXP(X): real X : Returns the natural logarithm of X.